classdef TriangleQuadrature < Quadrature
    properties 
        quadpts
        weights
    end
    methods 
        function obj = TriangleQuadrature(order)
            switch order
                case 1     % Order 1, nQuad 1
                    obj.quadpts = [1/3, 1/3, 1/3];
                    obj.weighs = 1;
                case 2     % Order 2, nQuad 3
                    obj.quadpts = [2/3, 1/6, 1/6; ...
                              1/6, 2/3, 1/6; ...
                              1/6, 1/6, 2/3];
                    obj.weights = [1/3, 1/3, 1/3];
                case 3     % Order 3, nQuad 4
                    obj.quadpts = [1/3, 1/3, 1/3; ...
                              0.6, 0.2, 0.2; ...
                              0.2, 0.6, 0.2; ...
                              0.2, 0.2, 0.6];
                    obj.weights = [-27/48, 25/48, 25/48, 25/48];
                case 4     % Order 4, nQuad 6
                    obj.quadpts = [0.108103018168070, 0.445948490915965, 0.445948490915965; ...
                              0.445948490915965, 0.108103018168070, 0.445948490915965; ...
                              0.445948490915965, 0.445948490915965, 0.108103018168070; ...
                              0.816847572980459, 0.091576213509771, 0.091576213509771; ...
                              0.091576213509771, 0.816847572980459, 0.091576213509771; ...
                              0.091576213509771, 0.091576213509771, 0.816847572980459];
                     obj.weights = [0.223381589678011, 0.223381589678011, 0.223381589678011, ...
                              0.109951743655322, 0.109951743655322, 0.109951743655322];
                case 5     % Order 5, nQuad 7
                    alpha1 = 0.059715871789770;      beta1 = 0.470142064105115;
                    alpha2 = 0.797426985353087;      beta2 = 0.101286507323456;
                    obj.quadpts = [   1/3,    1/3,    1/3; ...
                              alpha1,  beta1,  beta1; ...
                               beta1, alpha1,  beta1; ...
                               beta1,  beta1, alpha1; ...
                              alpha2,  beta2,  beta2; ...
                               beta2, alpha2,  beta2; ...
                               beta2,  beta2, alpha2];
                    obj.weights= [0.225, 0.132394152788506, 0.132394152788506, 0.132394152788506, ...
                        0.125939180544827, 0.125939180544827, 0.125939180544827];
                otherwise
                    error('I don not have such a quadrature!')
            end
        end
    end
end
